---
output: 
  pdf_document:
    citation_package: natbib
    keep_tex: true
    fig_caption: true
    latex_engine: pdflatex
    template: ~/Dropbox/Software/templates/svm-r-markdown-templates-master/svm-latex-ms.tex
# for template details, see http://svmiller.com/blog/2016/02/svm-r-markdown-manuscript/    
title: "Decomposing Public Opinion Variation into Ideology, Idiosyncrasy and Instability"
thanks: "**This version**: `r format(Sys.time(), '%B %d, %Y')`, forthcoming, Journal of Politics. Replication materials at doi:10.7910/DVN/KHBDWU"  
author:
- name: Benjamin E Lauderdale
  affiliation: London School of Economics and Political Science
- name: Chris Hanretty
  affiliation: Royal Holloway, University of London
- name: Nick Vivyan
  affiliation: Durham University
abstract: "We propose a method for decomposing variation in the issue preferences that US citizens express on surveys into three sources of variability that correspond to major threads in public opinion research.  We find that, averaging across a set of high profile US political issues, a single ideological dimension accounts for about 1/7 of opinion variation, individuals' idiosyncratic preferences account for about 3/7, and response instability for the remaining 3/7. These shares vary substantially across issue types and the average share attributable to ideology doubles when a second ideological dimension is permitted. We also find that (unidimensional) ideology accounts for almost twice as much response variation (and response instability is substantially lower) among respondents with high, rather than low, political knowledge.  Our estimation strategy is based on an ordinal probit model with random effects, and is applicable to other data sets that include repeated measurements of ordinal issue position data."
keywords: ""
date: "`r format(Sys.time(), '%B %d, %Y')`"
geometry: margin=1in
fontfamily: mathpazo
fontsize: 11pt
# papersize: a4paper
spacing: double
bibliography: master.bib
biblio-style: apsr
---

```{r, include=FALSE}

# Libraries

library(knitr)
library(xtable)
library(rstan)
library(abind)
library(ggplot2)
library(ggtern)
library(ggrepel)
library(MCMCpack)
library(reshape2)
library(dplyr)


# Settings

set.seed(42)
REFIT <- FALSE # Set to TRUE to fit models, FALSE to generate report from saved models
HMCSample <- 2000
HMCBurn <- 500
HMCChains <- 2
rstan_options(auto_write = TRUE)
options(mc.cores = HMCChains)

# Function definitions

rmse <- function(x1,x2) sqrt(mean((x1-x2)^2))
extract_chain <- function(var) extract(posterior,var)[[1]]
fadecol <- function(colours,alpha=1){
  rgb.dat <- col2rgb(colours)
  return(rgb(rgb.dat[1,],rgb.dat[2,],rgb.dat[3,],alpha*255,maxColorValue=255))
}  

```

```{r, include=FALSE}

# load survey data, see doi:10.7910/DVN/HOH1F0 for original data package and documentation
alldata <- read.csv("ssi-spatial-all-merged-with-wave1.tab",sep="\t")
nobs_full <- nrow(alldata)

# Select only respondents interviewed twice...

panelobs <- rowSums(!is.na(alldata[,c("eq_marijuana_POST","eq_environment_POST","eq_socialsecurity_POST",
                 "eq_guns_POST","eq_health_POST","eq_immigration_POST","eq_taxes_POST",
                 "eq_abortion_POST","eq_medicare_POST","eq_gays_POST","eq_unions_POST",
                 "eq_contraception_POST","eq_education_POST")])) != 0

t.test(alldata$pkscale ~ panelobs)
alldata <- subset(alldata,subset=panelobs)
nobs_panel <- nrow(alldata)

alldata$pkscale2 <- (alldata$pkscale - min(alldata$pkscale))/(max(alldata$pkscale) - min(alldata$pkscale))
varnames <- c("Marijuana","Environment","Social_Security","Guns","Health","Immigration","Taxes","Abortion","Medicare","Gay_Rights","Unions","Contraception","Education")
alldata$partycol <- c("#0000FF","#0000DD","#0000BB","#BB00BB","#BB0000","#DD0000","#FF0000")[alldata$partyid]

# Stan Model Data
Y1 <- alldata[,c("eq_marijuana","eq_environment","eq_socialsecurity",
                "eq_guns","eq_health","eq_immigration","eq_taxes",
                "eq_abortion","eq_medicare","eq_gays","eq_unions",
                "eq_contraception","eq_education")]
for (j in 1:ncol(Y1)) Y1[,j] <- as.integer(Y1[,j])
Y1 <- as.matrix(Y1)
colnames(Y1) <- varnames
Y1 <- replace(Y1,is.na(Y1),-1L)

Y2 <- alldata[,c("eq_marijuana_POST","eq_environment_POST","eq_socialsecurity_POST",
                 "eq_guns_POST","eq_health_POST","eq_immigration_POST","eq_taxes_POST",
                 "eq_abortion_POST","eq_medicare_POST","eq_gays_POST","eq_unions_POST",
                 "eq_contraception_POST","eq_education_POST")]
for (j in 1:ncol(Y2)) Y2[,j] <- as.integer(Y2[,j])
Y2 <- as.matrix(Y2)
colnames(Y2) <- varnames
Y2 <- replace(Y2,is.na(Y2),-1L)

Y <- abind(Y1,Y2,along=3)



# load legislator data
legislatordata <- read.csv("state-legislator-survey.tab",sep="\t")
Yleg <- legislatordata[,c("eq_environment",
                "eq_guns","eq_health","eq_immigration","eq_taxes",
                "eq_abortion","eq_medicare","eq_gays","eq_unions",
                "eq_education")]
Yleg <- replace(Yleg,is.na(Yleg),-1L)
# make comparable survey data from wave 1
Ylegcomp <- alldata[,c("eq_environment",
                "eq_guns","eq_health","eq_immigration","eq_taxes",
                "eq_abortion","eq_medicare","eq_gays","eq_unions",
                "eq_education")]
Ylegcomp <- replace(Ylegcomp,is.na(Ylegcomp),-1L)

```

# Introduction

Since @converse1964nature political scientists have debated the extent
to which citizens have organized preferences about political issues,
or indeed whether they have such preferences at all [@achen1975mass;
@zaller1992simple; @zaller1992nature].  In this short paper, we
exploit recent survey data collected by @broockman2016approaches to decompose
variation in citizens' expressed issue preferences into three sources
of variability.  First, how much is *ideological*, by which we mean
predictable from a citizen's preferences on other issues given the way
preferences are typically organized in the population?[^constraint]  Second, how
much is *idiosyncratic*, by which we mean preferences that are not
predictable from a citizen's preferences on other issues, but which
are nonetheless stable in repeated measurement?  Third, how much is
*instability*, by which we mean preference variation that is not
stable in repeated measurement and which combines measurement error
and non-attitudinal response behaviours?  Although scholars have at
various points argued for the relative importance of each of these
sources of variability in opinion, to our knowledge none have provided
a quantitative summary of how much variation is attributable to each
one.

[^constraint]: This is a thin definition of *ideology*---what Converse refers to as *constraint*. 

# Context

There is a long-running debate within political science about the
extent to which most citizens have "reasonably well-formed attitudes
on major political issues" [@zaller1992simple, 579] which exist
independently of the survey instruments designed to measure such
attitudes.  On one view, survey respondents may be perfectly capable
of responding to questions designed to elicit political preferences,
but these responses merely reflect the result of a sample of
considerations bearing on the matter, which are context-biased in ways
that lead them to give different answers at different points in time
[@zaller1992simple, 579]. Evidence in favour of this view comes from
the low correlation between many measures of attitudes at different
points in time; from open-ended survey questions eliciting voter
considerations; and from the extensive literature on priming.  Against
this evidence, it has been argued that low over-time correlation is an
indication that survey responses are subject to large measurement
error such that stable underlying attitudes can result in volatile
survey responses [@achen1975mass], and that many priming effects are
better understood as information effects [@lenz2009priming].

If citizens do not have stable preferences on major political issues
at all, then there is no point arguing about whether those attitudes
can be effectively summarized by a position in some latent ideological
space. However, if political preferences are at least partially stable
or real, then they may or may not have a low dimensional structure
(ideology) that explains much of the variation in preferences across
different issues. There is disagreement about whether preferences can
be represented using just one left-right dimension
[@jacoby1995structure; @jost2009political] or whether they must
instead be represented using two dimensions [@duckitt2009dual;
@feldman2014understanding; @evans1996measuring]. The more that
preferences can be represented in terms of a small number of
underlying dimensions, the less role there is for idiosyncratic
preferences, understood as `non-ideological' preferences which cannot
be predicted on the basis of preferences about other issues, but which
are nonetheless real, stable views about particular issues
[@broockman2016approaches].

This paper argues that we can usefully understand these two
debates---about the extent to which preferences are real and about the
extent to which they are ideological---as part of a larger question
about the composition of survey response variation. Our decomposition
of preferences into components corresponding to ideology, idiosyncrasy
and instability can be mapped onto these previous debates.  The debate
over whether preferences are real is a debate about the relative
contribution of ideology plus idiosyncrasy versus instability.  The
debate over whether preferences are ideological is a debate about the
relative contribution of ideology versus idiosyncrasy (either with or
without instability).

In our reading, no one has directly assessed the three-way
decomposition quantitatively.  Converse evaluates both the correlation
structure across issues as well as the stability of responses within
issues, but the synthesis of the evidence from those analyses is
qualitative [@converse1964nature].  Subsequent studies, including the
Broockman study whose data we re-analyse, have similarly done a series
of analyses and then summarized the results qualitatively.  The costs
of a purely qualitative synthesis of the evidence are apparent in the
subsequent multi-decade debate spurred by Converse's study, a debate
that tended to get simplified into binary questions about whether
preferences were real or not, or ideological or not, rather than
decompositional questions about how much of observed preference
variation is real and how much is ideological [@converse2000review;
@feldman2013political].


# Data

To quantitatively decompose response variation into ideology,
idiosyncrasy and instability, we need data with several
features. Repeated measurements on the same respondents must be taken
in order to distinguish stable preferences from response
instability. Questions on a range of political issues must be asked in
order to assess whether preferences are ideologically structured or
not.  Anchored response scales are desirable in order to limit
measurement error associated with differential scale use, both across
respondents and within respondents across repeated measurements.


Fortunately, data satisfying all of these criteria are provided by
@broockman2016approaches. Broockman surveyed US citizens in early
January and again in late February 2014. He asked eight questions as
part of a political knowledge battery and thirteen ordinal questions
on particular policies. The questions covered both economic issues
(health care, taxes, Medicare, unions, education, social security) and
social/cultural issues (gun control, immigration, abortion, the
environment, gay rights, affirmative action, and contraception). The
responses to the ordinal policy questions were constructed so that the
middle option represented the status quo, and the next least extreme
response options captured the median position of the Democratic and
Republican parties respectively. Remaining options represented more
extreme positions.[^responseorder] Giving such specific response
options mitigates the risk that variation in manifest survey responses
will result from the way respondents use response scales rather than
from respondents' varying underlying preferences.

[^responseorder]: Response options were ordered from left- to
    right-wing for some questions, and right- to left-wing for others.

```{r, eval = FALSE, echo = FALSE, results = "asis"}
box1 <- read.csv("response_example_3.csv",
                 fileEncoding = "utf-8")
print(xtable(box1[,-1,drop = FALSE],
             caption = "Example of a seven-point response scale for the prompt: Which statement comes closest to describing your views on taxes?",
             align = c("r","p{14cm}")),
      include.rownames = TRUE,
      booktabs = FALSE,
      comment=FALSE, floating=TRUE,table.placement="t")
```

Though we focus on the policy questions, we also use Broockman's
latent trait measure of political knowledge. Correct or incorrect
answers to eight knowledge questions were aggregated allowing for
variation in question difficulty and the ability of each question to
discriminate between respondents with otherwise similar levels of
knowledge. For simplicity, we categorise respondents as 'high' or
'low' political knowledge using the median score political knowledge
score as the cutoff. `r nobs_full` respondents participated in the first wave, of whom 
`r nobs_panel` also participated in the second wave.  Average levels of
political knowledge were were modestly (0.15 SD) higher among
wave two respondents than among respondents who only responded to wave
one.


# Methods

Decomposing variation in ordinal response data requires us to make
some scaling and modelling assumptions.  We assume an ordinal probit
response model, where respondents have latent continuous responses
($Y^*$) which map on to manifest categorical responses ($Y$) according to
cutpoints, which are drawn from a uniform
prior subject to the ordering constraint.  Our model for the latent
response takes the form:
$$
Y^*_{ijt} = \beta_j \theta_i + \nu_{ij} + \epsilon_{ijt}
$$
where $i$ indexes respondents, where $j$ indexes questions, where $t$
indexes the survey wave ($t = 1, 2$) and where the parameters in the
model are drawn $\theta_i \sim N(0, 1)$, $\nu_{ij} \sim N(0,
\omega_j^2)$, and $\epsilon_{ijt} \sim N(0, \sigma_j^2)$.

The first term in the model ($\beta_j\theta_i$) is the *ideological*
component, operationalised as a unidimensional spatial model where
$\beta_j$ describes how the respondent's latent ideal point $\theta$
predicts her latent response on issue $j$. The signs of $\theta$ and
$\beta$ are both arbitrary and unimportant for the variance
decomposition. The second term in the model ($\nu_{ij}$) is a
respondent-by-issue random effect. This term captures *idiosyncratic*
variation in preferences which is not attributable to ideology, but
which is nonetheless stable across repeat measurements in different time periods $t$. The third
term in the model ($\epsilon_{ijt}$) is an error term capturing
*instability*, variation in an individuals' manifest responses across repeated
measurements.  This *instability* could include response changes from several hypothesized mechanisms, including instability as sheer error [@achen1975mass], as sampling from top-of-the-head considerations [@zaller1992simple], as "flipping a coin" [@converse1964nature], or as `real' opinion change in the brief time between survey waves (which we discuss further in the conclusion).

We are not interested in the sign or magnitude of these terms for any
individual respondent, but rather in their variance per survey-item
across respondents. In the case of the respondent-by-issue random
effect and the error term, we use $\omega_j^2$ and $\sigma_j^2$ to
refer to these variances. The variance of the ideological component is
$\beta_j^2$, which follows from the identity $Var(cx) = c^2Var(x)$ and
the fact that $\theta_i \sim N(0, 1)$.

Because these sources of variation are independent by assumption, the
total variance per survey-item can be expressed as $\Sigma_j =
\beta_j^2 + \omega_j^2 + \sigma_j^2$. The share of the total variance
attributable to ideology can therefore be expressed as
$\beta_j^2/\Sigma_j$; the share attributable to idiosyncrasy as
$\omega_j^2/\Sigma_j$, and that attributable to instability
$\sigma_j^2/\Sigma_j$. This three-way decomposition of latent scale
variance is similar in spirit to the @mckelvey1975statistical
pseudo-R$^2$ statistic for the ordinal probit regression model, which
estimates the R$^2$ that would be recovered via linear regression on
the unobserved latent variable.

The priors for $\beta_j$, $\omega_j$ and $\sigma_j$ are standard
normal, half-normal, and half-normal, respectively.  These priors are
chosen so that the prior over the decomposition
$(\beta_j^2/\Sigma_j,\omega_j^2/\Sigma_j,\sigma_j^2/\Sigma_j)$ for
each issue $j$ is a symmetric Dirichlet distribution with parameters
equal to $1/2$, which is Jeffrey's non-informative
prior.[^priorchoice] These priors correspond to a weak prior
expectation for an equal split of the variance between ideology,
idiosyncrasy and instability, and any deviation from this in our results reflects evidence from the data. The model was estimated using Stan 2.16.2 [@stan].

[^priorchoice]: If $\beta_j$, $\omega_j$ and $\sigma_j$ have prior distributions
    that are standard (half-)normal, their squares are distributed
    $\chi_1^2$, which is equivalent to $\Gamma(1/2,2)$.  The simplex
    constructed by dividing three $\Gamma(1/2,2)$ random variables by
    their sum is distributed $D(1/2,1/2,1/2)$.
    
```{r, eval=FALSE, echo=FALSE}

# If you want to convince yourself of this by simulation...

draws <- 10000
x <- matrix(rnorm(3*draws),draws,3)
priordraws <- sweep(x^2,1,rowSums(x^2),"/")

dirchdraws <- rdirichlet(draws, c(0.5,0.5,0.5))

par(mfrow=c(1,2))
plot(decomp[,1],decomp[,2],col=rgb(0,0,0,0.1),pch=16,cex=0.5,main="Decomposition of Squares\nof Three Standard Normals",xlab="First Component",ylab="Second Component")
plot(dirchdraws[,1],dirchdraws[,2],col=rgb(0,0,0,0.1),pch=16,cex=0.5,main="Dirichlet (1/2,1/2,1/2)",xlab="First Component",ylab="Second Component")

```
	


```{r, include=FALSE}

standata <- list(
  Y = Y, # Responses to 7 pt scales
  X = as.integer(cut(alldata$pkscale2,c(-Inf,median(alldata$pkscale2),1))) ,
  N_respondents = nrow(Y),
  N_items = ncol(Y),
  N_categories = max(Y)
)


standataleg <- list(
  Y = as.matrix(Yleg), # Responses to 7 pt scales
  N_respondents = nrow(Yleg),
  N_items = ncol(Yleg),
  N_categories = max(Yleg)
)


standatalegcomp <- list(
  Y = as.matrix(Ylegcomp), # Responses to 7 pt scales
  N_respondents = nrow(Ylegcomp),
  N_items = ncol(Ylegcomp),
  N_categories = max(Ylegcomp)
)

# Stan Model Code

## Model without spatial dimension

stancode0 <- '

data {
  int<lower = 1> N_respondents; // Number of respondents 
  int<lower = 1> N_items; // Number of questions
  int<lower = 1> N_categories; // Number of response categories
  int<lower = -1,upper=N_categories> Y[N_respondents,N_items,2]; // Array of response data  
}

parameters {
  vector<lower=0>[N_items] sigma_std;
  vector<lower=0>[N_items] omega_std;
  ordered[N_categories-1] alpha[N_items];
  real nu_std[N_respondents,N_items]; // person-item random effects
}

transformed parameters {

  real<lower=0> norm[N_items];
  real<lower=0> omega[N_items];
  real<lower=0> sigma[N_items];
  vector[N_items] beta;
  real nu[N_respondents,N_items]; // person-item random effects
  real eta[N_respondents,N_items];
  vector[N_categories] pi[N_respondents,N_items];

  for(j in 1:N_items) {  // Normalize so that latent scale has variance equal to one for each item
    norm[j] <- sigma_std[j]*sigma_std[j] + omega_std[j]*omega_std[j];
    sigma[j] <- sigma_std[j]/sqrt(norm[j]); 
    omega[j] <- omega_std[j]/sqrt(norm[j]); 
    beta[j] <- 0.0; 
  }

  for(i in 1:N_respondents){
    for(j in 1:N_items){
      nu[i,j] <- nu_std[i,j]*omega[j];
      eta[i,j] <- nu[i,j]; 
      pi[i,j,1] <- 1.0 - Phi((eta[i,j] - alpha[j,1])/sigma[j]);
      for (k in 2:(N_categories-1)) pi[i,j,k] <- Phi((eta[i,j] - alpha[j,k-1])/sigma[j]) - Phi((eta[i,j] - alpha[j,k])/sigma[j]);
      pi[i,j,N_categories] <- Phi((eta[i,j] - alpha[j,N_categories-1])/sigma[j]);
    }
  }
}

model {

  for(j in 1:N_items){
    for(i in 1:N_respondents){
      nu_std[i,j] ~ normal(0,1);
      if (Y[i,j,1] != -1) Y[i,j,1] ~ categorical(pi[i,j]);
      if (Y[i,j,2] != -1) Y[i,j,2] ~ categorical(pi[i,j]);
    }
  }
  omega_std ~ normal(0,1);
  sigma_std ~ normal(0,1);
}

generated quantities {

  simplex[3] decomposition[N_items];
  for(j in 1:N_items) {
    decomposition[j,1] <- beta[j]*beta[j];
    decomposition[j,2] <- omega[j]*omega[j];
    decomposition[j,3] <- sigma[j]*sigma[j];
  }

}
'


## Model without second wave (and thus instability estimate)

stancode0b <- '

data {
  int<lower = 1> N_respondents; // Number of respondents 
  int<lower = 1> N_items; // Number of questions
  int<lower = 1> N_categories; // Number of response categories
  int<lower = -1,upper=N_categories> Y[N_respondents,N_items]; // Array of response data  
}

parameters {
  vector[N_items] beta_std;
  vector<lower=0>[N_items] omega_std;
  vector<lower=0>[N_items] sigma_std;
  vector[N_respondents] theta;
  ordered[N_categories-1] alpha[N_items];
  real nu_std[N_respondents,N_items]; // person-item random effects

}

transformed parameters {

  real<lower=0> norm[N_items];
  real<lower=0> omega[N_items];
  real<lower=0> sigma[N_items];
  vector[N_items] beta;
  real nu[N_respondents,N_items]; // person-item random effects
  real eta[N_respondents,N_items];
  vector[N_categories] pi[N_respondents,N_items];

  for(j in 1:N_items) {  // Normalize so that latent scale has variance equal to one for each item
    norm[j] <- sigma_std[j]*sigma_std[j] + omega_std[j]*omega_std[j] + beta_std[j]*beta_std[j];
    sigma[j] <- sigma_std[j]/sqrt(norm[j]); 
    omega[j] <- omega_std[j]/sqrt(norm[j]); 
    beta[j] <- beta_std[j]/sqrt(norm[j]);
  }

  for(i in 1:N_respondents){
    for(j in 1:N_items){
      nu[i,j] <- nu_std[i,j]*omega[j];
      eta[i,j] <- beta[j]*theta[i] + nu[i,j]; // spatial model
      pi[i,j,1] <- 1.0 - Phi((eta[i,j] - alpha[j,1])/sigma[j]);
      for (k in 2:(N_categories-1)) pi[i,j,k] <- Phi((eta[i,j] - alpha[j,k-1])/sigma[j]) - Phi((eta[i,j] - alpha[j,k])/sigma[j]);
      pi[i,j,N_categories] <- Phi((eta[i,j] - alpha[j,N_categories-1])/sigma[j]);
    }
  }
}

model {

  for(j in 1:N_items){
    for(i in 1:N_respondents){
      nu_std[i,j] ~ normal(0,1);
      if (Y[i,j] != -1) Y[i,j] ~ categorical(pi[i,j]);
    }
  }
  omega_std ~ normal(0,1);
  sigma_std ~ normal(0,1);
  beta_std ~ normal(0,1);
  theta ~ normal(0,1);
}

generated quantities {

  simplex[3] decomposition[N_items];
  for(j in 1:N_items) {
    decomposition[j,1] <- beta[j]*beta[j];
    decomposition[j,2] <- omega[j]*omega[j];
    decomposition[j,3] <- sigma[j]*sigma[j];
  }

}

'



## Model with spatial dimension

stancode1 <- '

data {
  int<lower = 1> N_respondents; // Number of respondents 
  int<lower = 1> N_items; // Number of questions
  int<lower = 1> N_categories; // Number of response categories
  int<lower = -1,upper=N_categories> Y[N_respondents,N_items,2]; // Array of response data  
}

parameters {
  vector[N_items] beta_std;
  vector<lower=0>[N_items] omega_std;
  vector<lower=0>[N_items] sigma_std;
  vector[N_respondents] theta;
  ordered[N_categories-1] alpha[N_items];
  real nu_std[N_respondents,N_items]; // person-item random effects

}

transformed parameters {

  real<lower=0> norm[N_items];
  real<lower=0> omega[N_items];
  real<lower=0> sigma[N_items];
  vector[N_items] beta;
  real nu[N_respondents,N_items]; // person-item random effects
  real eta[N_respondents,N_items];
  vector[N_categories] pi[N_respondents,N_items];

  for(j in 1:N_items) {  // Normalize so that latent scale has variance equal to one for each item
    norm[j] <- sigma_std[j]*sigma_std[j] + omega_std[j]*omega_std[j] + beta_std[j]*beta_std[j];
    sigma[j] <- sigma_std[j]/sqrt(norm[j]); 
    omega[j] <- omega_std[j]/sqrt(norm[j]); 
    beta[j] <- beta_std[j]/sqrt(norm[j]);
  }

  for(i in 1:N_respondents){
    for(j in 1:N_items){
      nu[i,j] <- nu_std[i,j]*omega[j];
      eta[i,j] <- beta[j]*theta[i] + nu[i,j]; // spatial model
      pi[i,j,1] <- 1.0 - Phi((eta[i,j] - alpha[j,1])/sigma[j]);
      for (k in 2:(N_categories-1)) pi[i,j,k] <- Phi((eta[i,j] - alpha[j,k-1])/sigma[j]) - Phi((eta[i,j] - alpha[j,k])/sigma[j]);
      pi[i,j,N_categories] <- Phi((eta[i,j] - alpha[j,N_categories-1])/sigma[j]);
    }
  }
}

model {

  for(j in 1:N_items){
    for(i in 1:N_respondents){
      nu_std[i,j] ~ normal(0,1);
      if (Y[i,j,1] != -1) Y[i,j,1] ~ categorical(pi[i,j]);
      if (Y[i,j,2] != -1) Y[i,j,2] ~ categorical(pi[i,j]);
    }
  }
  omega_std ~ normal(0,1);
  sigma_std ~ normal(0,1);
  beta_std ~ normal(0,1);
  theta ~ normal(0,1);
}

generated quantities {

  simplex[3] decomposition[N_items];
  for(j in 1:N_items) {
    decomposition[j,1] <- beta[j]*beta[j];
    decomposition[j,2] <- omega[j]*omega[j];
    decomposition[j,3] <- sigma[j]*sigma[j];
  }

}

'

## Model with two spatial dimensions

stancode1b <- '

data {
  int<lower = 1> N_respondents; // Number of respondents 
  int<lower = 1> N_items; // Number of questions
  int<lower = 1> N_categories; // Number of response categories
  int<lower = -1,upper=N_categories> Y[N_respondents,N_items,2]; // Array of response data  
}

parameters {
  vector[N_items] beta_std[2];
  vector<lower=0>[N_items] omega_std;
  vector<lower=0>[N_items] sigma_std;
  vector[N_respondents] theta[2];
  ordered[N_categories-1] alpha[N_items];
  real nu_std[N_respondents,N_items]; // person-item random effects

}

transformed parameters {

  real<lower=0> norm[N_items];
  real<lower=0> omega[N_items];
  real<lower=0> sigma[N_items];
  vector[N_items] beta[2];
  real nu[N_respondents,N_items]; // person-item random effects
  real eta[N_respondents,N_items];
  vector[N_categories] pi[N_respondents,N_items];

  for(j in 1:N_items) {  // Normalize so that latent scale has variance equal to one for each item
    norm[j] <- sigma_std[j]*sigma_std[j] + omega_std[j]*omega_std[j] + beta_std[1,j]*beta_std[1,j] + beta_std[2,j]*beta_std[2,j];
    sigma[j] <- sigma_std[j]/sqrt(norm[j]); 
    omega[j] <- omega_std[j]/sqrt(norm[j]); 
    beta[1,j] <- beta_std[1,j]/sqrt(norm[j]); 
    beta[2,j] <- beta_std[2,j]/sqrt(norm[j]); 
  }

  for(i in 1:N_respondents){
    for(j in 1:N_items){
      nu[i,j] <- nu_std[i,j]*omega[j];
      eta[i,j] <- beta[1,j]*theta[1,i] + beta[2,j]*theta[2,i] + nu[i,j]; // spatial model
      pi[i,j,1] <- 1.0 - Phi((eta[i,j] - alpha[j,1])/sigma[j]);
      for (k in 2:(N_categories-1)) pi[i,j,k] <- Phi((eta[i,j] - alpha[j,k-1])/sigma[j]) - Phi((eta[i,j] - alpha[j,k])/sigma[j]);
      pi[i,j,N_categories] <- Phi((eta[i,j] - alpha[j,N_categories-1])/sigma[j]);
    }
  }
}

model {

  for(j in 1:N_items){
    for(i in 1:N_respondents){
      nu_std[i,j] ~ normal(0,1);
      if (Y[i,j,1] != -1) Y[i,j,1] ~ categorical(pi[i,j]);
      if (Y[i,j,2] != -1) Y[i,j,2] ~ categorical(pi[i,j]);
    }
  }
  omega_std ~ normal(0,1);
  sigma_std ~ normal(0,1);
  beta_std[1] ~ normal(0,sqrt(1.0/2.0));
  beta_std[2] ~ normal(0,sqrt(1.0/2.0));
  theta[1] ~ normal(0,1);
  theta[2] ~ normal(0,1);
}

generated quantities {

  simplex[3] decomposition[N_items];
  for(j in 1:N_items) {
    decomposition[j,1] <- beta[1,j]*beta[1,j]+beta[2,j]*beta[2,j];
    decomposition[j,2] <- omega[j]*omega[j];
    decomposition[j,3] <- sigma[j]*sigma[j];
  }

}

'



## Model with three spatial dimensions (Fuck everything, we are doing five blades.)

stancode1c <- '

data {
  int<lower = 1> N_respondents; // Number of respondents 
  int<lower = 1> N_items; // Number of questions
  int<lower = 1> N_categories; // Number of response categories
  int<lower = -1,upper=N_categories> Y[N_respondents,N_items,2]; // Array of response data  
}

parameters {
  vector[N_items] beta_std[3];
  vector<lower=0>[N_items] omega_std;
  vector<lower=0>[N_items] sigma_std;
  vector[N_respondents] theta[3];
  ordered[N_categories-1] alpha[N_items];
  real nu_std[N_respondents,N_items]; // person-item random effects

}

transformed parameters {

  real<lower=0> norm[N_items];
  real<lower=0> omega[N_items];
  real<lower=0> sigma[N_items];
  vector[N_items] beta[3];
  real nu[N_respondents,N_items]; // person-item random effects
  real eta[N_respondents,N_items];
  vector[N_categories] pi[N_respondents,N_items];

  for(j in 1:N_items) {  // Normalize so that latent scale has variance equal to one for each item
    norm[j] <- sigma_std[j]*sigma_std[j] + omega_std[j]*omega_std[j] + beta_std[1,j]*beta_std[1,j] + beta_std[2,j]*beta_std[2,j] + beta_std[3,j]*beta_std[3,j];
    sigma[j] <- sigma_std[j]/sqrt(norm[j]); 
    omega[j] <- omega_std[j]/sqrt(norm[j]); 
    beta[1,j] <- beta_std[1,j]/sqrt(norm[j]); 
    beta[2,j] <- beta_std[2,j]/sqrt(norm[j]); 
    beta[3,j] <- beta_std[3,j]/sqrt(norm[j]); 
  }

  for(i in 1:N_respondents){
    for(j in 1:N_items){
      nu[i,j] <- nu_std[i,j]*omega[j];
      eta[i,j] <- beta[1,j]*theta[1,i] + beta[2,j]*theta[2,i] + beta[3,j]*theta[3,i] + nu[i,j]; // spatial model
      pi[i,j,1] <- 1.0 - Phi((eta[i,j] - alpha[j,1])/sigma[j]);
      for (k in 2:(N_categories-1)) pi[i,j,k] <- Phi((eta[i,j] - alpha[j,k-1])/sigma[j]) - Phi((eta[i,j] - alpha[j,k])/sigma[j]);
      pi[i,j,N_categories] <- Phi((eta[i,j] - alpha[j,N_categories-1])/sigma[j]);
    }
  }
}

model {

  for(j in 1:N_items){
    for(i in 1:N_respondents){
      nu_std[i,j] ~ normal(0,1);
      if (Y[i,j,1] != -1) Y[i,j,1] ~ categorical(pi[i,j]);
      if (Y[i,j,2] != -1) Y[i,j,2] ~ categorical(pi[i,j]);
    }
  }
  omega_std ~ normal(0,1);
  sigma_std ~ normal(0,1);
  beta_std[1] ~ normal(0,sqrt(1.0/3.0));
  beta_std[2] ~ normal(0,sqrt(1.0/3.0));
  beta_std[3] ~ normal(0,sqrt(1.0/3.0));
  theta[1] ~ normal(0,1);
  theta[2] ~ normal(0,1);
  theta[3] ~ normal(0,1);
}

generated quantities {

  simplex[3] decomposition[N_items];
  for(j in 1:N_items) {
    decomposition[j,1] <- beta[1,j]*beta[1,j]+beta[2,j]*beta[2,j]+beta[3,j]*beta[3,j];
    decomposition[j,2] <- omega[j]*omega[j];
    decomposition[j,3] <- sigma[j]*sigma[j];
  }

}

'

## Model splitting by political knowledge

stancode2 <- '

data {
  int<lower = 1> N_respondents; // Number of respondents 
  int<lower = 1> N_items; // Number of questions
  int<lower = 1> N_categories; // Number of response categories
  int<lower = -1,upper=N_categories> Y[N_respondents,N_items,2]; // Array of response data  
  int<lower = 1,upper=2> X[N_respondents]; // Political mnowledge measure
}

parameters {
  real<lower=0> omega_std[N_items,2];
  real<lower=0> sigma_std[N_items,2];
  real beta_std[N_items,2];
  vector[N_respondents] theta;
  ordered[N_categories-1] alpha[N_items];
  real nu_std[N_respondents,N_items]; // person-item random effects
}

transformed parameters {

  real<lower=0> norm[N_items,2];
  real<lower=0> omega[N_items,2];
  real<lower=0> sigma[N_items,2];
  real beta[N_items,2];
  real nu[N_respondents,N_items]; // person-item random effects
  real eta[N_respondents,N_items];
  vector[N_categories] pi[N_respondents,N_items];

  for(k in 1:2){
    for(j in 1:N_items) {  // Normalize so that latent scale has variance equal to one for each item
      norm[j,k] <- sigma_std[j,k]*sigma_std[j,k] + omega_std[j,k]*omega_std[j,k] + beta_std[j,k]*beta_std[j,k];
      sigma[j,k] <- sigma_std[j,k]/sqrt(norm[j,k]); 
      omega[j,k] <- omega_std[j,k]/sqrt(norm[j,k]); 
      beta[j,k] <- beta_std[j,k]/sqrt(norm[j,k]); 
    }
  }

  for(i in 1:N_respondents){
    for(j in 1:N_items){
      nu[i,j] <- nu_std[i,j]*omega[j,X[i]];
      eta[i,j] <- beta[j,X[i]]*theta[i] + nu[i,j]; // spatial model
      pi[i,j,1] <- 1.0 - Phi((eta[i,j] - alpha[j,1])/sigma[j,X[i]]);
      for (k in 2:(N_categories-1)) pi[i,j,k] <- Phi((eta[i,j] - alpha[j,k-1])/sigma[j,X[i]]) - Phi((eta[i,j] - alpha[j,k])/sigma[j,X[i]]);
      pi[i,j,N_categories] <- Phi((eta[i,j] - alpha[j,N_categories-1])/sigma[j,X[i]]);
    }
  }
}

model {

  for(j in 1:N_items){
    for(i in 1:N_respondents){
      nu_std[i,j] ~ normal(0,1);
      if (Y[i,j,1] != -1) Y[i,j,1] ~ categorical(pi[i,j]);
      if (Y[i,j,2] != -1) Y[i,j,2] ~ categorical(pi[i,j]);
    }
  }
  for(j in 1:N_items){
    for (k in 1:2){
      omega_std[j,k] ~ normal(0,1);
      sigma_std[j,k] ~ normal(0,1);
      beta_std[j,k] ~ normal(0,1);
    }
  }
  theta ~ normal(0,1);
}

generated quantities {

  simplex[3] decomposition[N_items,2];
  for(j in 1:N_items) {
    decomposition[j,1,1] <- beta[j,1]*beta[j,1];
    decomposition[j,1,2] <- omega[j,1]*omega[j,1];
    decomposition[j,1,3] <- sigma[j,1]*sigma[j,1];
    decomposition[j,2,1] <- beta[j,2]*beta[j,2];
    decomposition[j,2,2] <- omega[j,2]*omega[j,2];
    decomposition[j,2,3] <- sigma[j,2]*sigma[j,2];
  }

}

'

geninit0 <- function() return(list(
  beta_std = rep(0,standata$N_items),
  sigma_std = rep(1,standata$N_items),
  omega_std = rep(1,standata$N_items),  
  theta = rnorm(standata$N_respondents,0,1),
  nu_std = matrix(rnorm(standata$N_respondents*standata$N_items,0,0.1),standata$N_respondents,standata$N_items),
  alpha = matrix(seq(-2,2,length.out=standata$N_categories-1),standata$N_items,standata$N_categories-1,byrow=TRUE)
))

geninit0bcit <- function() return(list(
  beta_std = rep(0,standatalegcomp$N_items),
  sigma_std = rep(1,standatalegcomp$N_items),
  omega_std = rep(1,standatalegcomp$N_items),  
  theta = rnorm(standatalegcomp$N_respondents,0,1),
  nu_std = matrix(rnorm(standatalegcomp$N_respondents*standatalegcomp$N_items,0,0.1),standatalegcomp$N_respondents,standatalegcomp$N_items),
  alpha = matrix(seq(-2,2,length.out=standatalegcomp$N_categories-1),standatalegcomp$N_items,standatalegcomp$N_categories-1,byrow=TRUE)
))

geninit0bleg <- function() return(list(
  beta_std = rep(0,standataleg$N_items),
  sigma_std = rep(1,standataleg$N_items),
  omega_std = rep(1,standataleg$N_items),  
  theta = rnorm(standataleg$N_respondents,0,1),
  nu_std = matrix(rnorm(standataleg$N_respondents*standataleg$N_items,0,0.1),standataleg$N_respondents,standataleg$N_items),
  alpha = matrix(seq(-2,2,length.out=standataleg$N_categories-1),standataleg$N_items,standataleg$N_categories-1,byrow=TRUE)
))

geninit1 <- function() return(list(
  beta_std = rep(0,standata$N_items),
  sigma_std = rep(1,standata$N_items),
  omega_std = rep(1,standata$N_items),  
  theta = rnorm(standata$N_respondents,0,1),
  nu_std = matrix(rnorm(standata$N_respondents*standata$N_items,0,0.1),standata$N_respondents,standata$N_items),
  alpha = matrix(seq(-2,2,length.out=standata$N_categories-1),standata$N_items,standata$N_categories-1,byrow=TRUE)
))

geninit1b <- function() return(list(
  beta_std = rbind(rep(0,standata$N_items),rep(0,standata$N_items)),
  sigma_std = rep(1,standata$N_items),
  omega_std = rep(1,standata$N_items),  
  theta = rbind(rnorm(standata$N_respondents,0,1),rnorm(standata$N_respondents,0,1)),
  nu_std = matrix(rnorm(standata$N_respondents*standata$N_items,0,0.1),standata$N_respondents,standata$N_items),
  alpha = matrix(seq(-2,2,length.out=standata$N_categories-1),standata$N_items,standata$N_categories-1,byrow=TRUE)
))

geninit1c <- function() return(list(
  beta_std = rbind(rep(0,standata$N_items),rep(0,standata$N_items),rep(0,standata$N_items)),
  sigma_std = rep(1,standata$N_items),
  omega_std = rep(1,standata$N_items),  
  theta = rbind(rnorm(standata$N_respondents,0,1),rnorm(standata$N_respondents,0,1),rnorm(standata$N_respondents,0,1)),
  nu_std = matrix(rnorm(standata$N_respondents*standata$N_items,0,0.1),standata$N_respondents,standata$N_items),
  alpha = matrix(seq(-2,2,length.out=standata$N_categories-1),standata$N_items,standata$N_categories-1,byrow=TRUE)
))

geninit2 <- function() return(list(
  beta_std = cbind(rep(0,standata$N_items),rep(0,standata$N_items)),
  sigma_std = cbind(rep(1,standata$N_items),rep(1,standata$N_items)),
  omega_std = cbind(rep(1,standata$N_items),rep(1,standata$N_items)),  
  theta = rnorm(standata$N_respondents,0,1),
  nu_std = matrix(rnorm(standata$N_respondents*standata$N_items,0,0.1),standata$N_respondents,standata$N_items),
  alpha = matrix(seq(-2,2,length.out=standata$N_categories-1),standata$N_items,standata$N_categories-1,byrow=TRUE)
))

# Estimate via Hamiltonian Monte Carlo

if (REFIT) {
  
  posterior0 <- stan(
    model_code = stancode0, 
    data = standata,
    init= geninit0,
    pars=c(
      "alpha","sigma","omega","decomposition"
    ), 
    iter = HMCSample + HMCBurn,
    warmup = HMCBurn, 
    chains = HMCChains,
    refresh = 1
  )  
  
  posterior0bcit <- stan(
    model_code = stancode0b, 
    data = standatalegcomp,
    init= geninit0bcit,
    pars=c(
      "alpha","sigma","omega","decomposition"
    ), 
    iter = HMCSample + HMCBurn,
    warmup = HMCBurn, 
    chains = HMCChains,
    refresh = 1
  )  

  posterior0bleg <- stan(
    model_code = stancode0b, 
    data = standataleg,
    init= geninit0bleg,
    pars=c(
      "alpha","sigma","omega","decomposition"
    ), 
    iter = HMCSample + HMCBurn,
    warmup = HMCBurn, 
    chains = HMCChains,
    refresh = 1
  )      
  
  posterior1 <- stan(
    model_code = stancode1, 
    data = standata,
    init= geninit1,
    pars=c(
      "alpha","beta","theta","sigma","omega","decomposition"
    ), 
    iter = HMCSample + HMCBurn,
    warmup = HMCBurn, 
    chains = HMCChains,
    refresh = 1
  )

  posterior1b <- stan(
    model_code = stancode1b, 
    data = standata,
    init= geninit1b,
    pars=c(
      "alpha","beta","theta","sigma","omega","decomposition"
    ), 
    iter = HMCSample + HMCBurn,
    warmup = HMCBurn, 
    chains = HMCChains,
    refresh = 1
  ) 
  
  posterior1c <- stan(
    model_code = stancode1c, 
    data = standata,
    init= geninit1c,
    pars=c(
      "alpha","beta","theta","sigma","omega","decomposition"
    ), 
    iter = HMCSample + HMCBurn,
    warmup = HMCBurn, 
    chains = HMCChains,
    refresh = 1
  )   
   
  posterior2 <- stan(
    model_code = stancode2, 
    data = standata,
    init= geninit2,
    pars=c(
      "alpha","beta","theta","sigma","omega","decomposition"
    ), 
    iter = HMCSample + HMCBurn,
    warmup = HMCBurn, 
    chains = HMCChains,
    refresh = 1
  )  
  
  save(posterior0,posterior0bcit,posterior0bleg,posterior1,posterior1b,posterior1c,posterior2,file = "stanout.Rdata")
  
} else load("stanout.Rdata")  

```


# Results

```{r, include=FALSE}

# Examine estimates

posterior <- posterior0
decomposition_chain_0 <- extract_chain("decomposition")
decomposition_est_0 <- t(apply(decomposition_chain_0,c(2,3),mean))
colnames(decomposition_est_0) <- varnames
rownames(decomposition_est_0) <- c("Ideology","Idiosyncrasy","Instability")
print(decomposition_est_0)


posterior <- posterior0bcit
decomposition_chain_0bcit <- extract_chain("decomposition")
decomposition_est_0bcit <- t(apply(decomposition_chain_0bcit,c(2,3),mean))
rownames(decomposition_est_0bcit) <- c("Ideology","Idiosyncrasy","Instability")
print(decomposition_est_0bcit)

posterior <- posterior0bleg
decomposition_chain_0bleg <- extract_chain("decomposition")
decomposition_est_0bleg <- t(apply(decomposition_chain_0bleg,c(2,3),mean))
rownames(decomposition_est_0bleg) <- c("Ideology","Idiosyncrasy","Instability")
print(decomposition_est_0bleg)

posterior <- posterior1
decomposition_chain_1 <- extract_chain("decomposition")
decomposition_est_1 <- t(apply(decomposition_chain_1,c(2,3),mean))
decomposition_se_1 <- t(apply(decomposition_chain_1,c(2,3),sd))
colnames(decomposition_est_1) <- colnames(decomposition_se_1) <- varnames
rownames(decomposition_est_1) <- rownames(decomposition_se_1) <- c("Ideology","Idiosyncrasy","Instability")
print(decomposition_est_1)

rms_difference_between_null_and_spatial_model_instability_shares <- 100*rmse(decomposition_est_0[3,],decomposition_est_1[3,])

posterior <- posterior1b
decomposition_chain_1b <- extract_chain("decomposition")
decomposition_est_1b <- t(apply(decomposition_chain_1b,c(2,3),mean))
decomposition_se_1b <- t(apply(decomposition_chain_1b,c(2,3),sd))
colnames(decomposition_est_1b) <- colnames(decomposition_se_1b) <- varnames
rownames(decomposition_est_1b) <- rownames(decomposition_se_1b) <- c("Ideology","Idiosyncrasy","Instability")
print(decomposition_est_1b)

posterior <- posterior1c
decomposition_chain_1c <- extract_chain("decomposition")
decomposition_est_1c <- t(apply(decomposition_chain_1c,c(2,3),mean))
decomposition_se_1c <- t(apply(decomposition_chain_1c,c(2,3),sd))
colnames(decomposition_est_1c) <- colnames(decomposition_se_1c) <- varnames
rownames(decomposition_est_1c) <- rownames(decomposition_se_1c) <- c("Ideology","Idiosyncrasy","Instability")
print(decomposition_est_1c)

posterior <- posterior2
decomposition_chain_2 <- extract_chain("decomposition")
decomposition_est_2 <- apply(decomposition_chain_2,c(2,3,4),mean)
decomposition_se_2 <- apply(decomposition_chain_2,c(2,3,4),sd)
dimnames(decomposition_est_2)[[1]] <- dimnames(decomposition_se_2)[[1]] <- varnames
dimnames(decomposition_est_2)[[2]] <- dimnames(decomposition_se_2)[[2]] <- c("Low Knowledge","High Knowledge")
dimnames(decomposition_est_2)[[3]] <- dimnames(decomposition_se_2)[[3]] <- c("Ideology","Idiosyncrasy","Instability")
print(decomposition_est_2)

overall_decomposition_table <- t(decomposition_est_1)*100
sortbyideo <- order(overall_decomposition_table[,1])
sortbyinst <- order(overall_decomposition_table[,3])
overall_decomposition_table <- rbind(overall_decomposition_table[sortbyideo,],colMeans(overall_decomposition_table)) 
rownames(overall_decomposition_table)[nrow(overall_decomposition_table)] <- "Average"
print(round(overall_decomposition_table[sortbyideo,]))

# Test differences between averages

decomposition_average_chain_1 <- apply(decomposition_chain_1,c(1,3),mean)
decomposition_average_int <- 100*apply(decomposition_average_chain_1,2,quantile,c(0.025,0.975))

decomposition_average_chain_1b <- apply(decomposition_chain_1b,c(1,3),mean)
decomposition_average_int_1b <- 100*apply(decomposition_average_chain_1b,2,quantile,c(0.025,0.975))

# Work out effect of dropping issues from average
drop_one_decomposition <- sapply(1:13, function(i) colMeans(apply(decomposition_chain_1[,-i,], c(1,3), mean)))
drop_one_ideo <- 100 * range(drop_one_decomposition[1,])

```

Figure 1 shows the estimated composition of opinion variation
for each issue in the data set that we examine. Averaging across 
issues, we estimate that the fraction of variation accounted for by
ideology is low 
(`r round(overall_decomposition_table[nrow(overall_decomposition_table),1])`%,
95% interval: 
`r floor(decomposition_average_int[1,1])`-`r ceiling(decomposition_average_int[2,1])`) when compared to that
explained by idiosyncrasy 
(`r round(overall_decomposition_table[nrow(overall_decomposition_table),2])`%,
95% interval: `r floor(decomposition_average_int[1,2])`-`r ceiling(decomposition_average_int[2,2])`) or instability 
(`r round(overall_decomposition_table[nrow(overall_decomposition_table),3])`%,
95% interval: `r floor(decomposition_average_int[1,3])`-`r ceiling(decomposition_average_int[2,3])`). Most of the stable
preference variation across individuals and issues is not organised
along a single political dimension.

The four most ideological issues are those related to the size of
government: Medicare, Taxes, Social Security and Health.  However, the
four issues with the most stable preferences are mostly different
issues: Immigration, Marijuana, Gay Rights and Health.  This
highlights the distinction between citizens having stable preferences
on issues and having ideological preferences.  Respondents to this
study have real and measurable views about marijuana and immigration
(relatively low instability) but those views are almost entirely
unpredictable on the basis of a unidimensional summary of their other
positions.  Stated preferences regarding unions (and to a
lesser extent education and the environment) are both lacking in
ideological structure and highly unstable; these are the issues on
which the most respondents seem to lack real preferences.

We have conducted a large number of robustness checks on variant forms of the model and data set (a full tabulation of these is provided in an online appendix).  When we simulate data from the estimates and data generating process of the model and refit the model, we recover the estimated parameter values and decomposition almost exactly, so there is negligible estimator bias.  If we omit the spatial term and decompose preference variation
into only idiosyncrasy and instability, we get nearly identical
estimates for instability, and the estimates of idiosyncrasy absorb the variation
associated with ideology in the full model. If we omit the idiosyncrasy term instead, the share estimated for ideology absorbs some of the stable variation, increasing from `r round(overall_decomposition_table[nrow(overall_decomposition_table),1])`% to 18%.

If we analyse only the first wave of data, we get nearly the same estimate for the share of
variation associated with ideology, but we can no longer distinguish
between idiosyncrasy and instability. For comparison purposes, this one-wave model can also be applied to Broockman's accompanying survey of state legislators. Averaging across the policy issues about which both state legislators and citizens were asked, we estimate that ideology accounts for `r round(rowMeans(decomposition_est_0bleg*100)[1])`% of the variation among state legislators, versus `r round(rowMeans(decomposition_est_0bcit*100)[1])`% for citizens.  This contrast illustrates that it is possible to estimate high shares of variation for ideology under this approach, as well as quantifying the well-known, very large difference between the degree of spatial structure in the positions of elites and citizens.

Adding a second spatial dimension to the ideological component of the model increases the fraction of variation explained by ideology from `r round(rowMeans(decomposition_est_1*100)[1])`% to `r round(rowMeans(decomposition_est_1b*100)[1])`% (95% interval: `r floor(decomposition_average_int_1b[1,1])`-`r ceiling(decomposition_average_int_1b[2,1])`) while leaving that due to instability unchanged. The fact that adding a second dimension nearly doubles the variance explained by ideology reinforces the point that a modest fraction of the variation in public opinion is captured by a single dimension, despite a two-party system that one might expect to organize opinion unidimensionally. 

Our estimates are somewhat sensitive to collapsing the extreme policy categories, with the ideology component rising to 16% when we collapse the two most extreme categories on each side to generate five point scales, and 20% when we collapse all categories on each side of the central status quo category to generate three point scales.  The increase in the ideology component is entirely at the expense of the inconsistency component, indicating that there is more structure and stability in the *direction* of change that respondents would like than in exactly *how far* they would like policy to move. Finally, we have analysed a version of the model where $\theta_i$ is allowed to change at the individual-level between the two periods as an alternative to the unstructured change captured by the instability term of the model, but we see very little evidence of ideologically coherent movement between these two survey waves.

```{r,results='asis',echo=FALSE, eval = FALSE}
print(xtable(overall_decomposition_table, digits=0, align=c("l","c","c","c"), caption="Percent of response variation attributable to ideology, idiosyncrasy and instability for each of 13 issues."), comment=FALSE, floating=TRUE,table.placement="t",hline.after = 13)
```


```{r, echo = FALSE, fig = TRUE, fig.cap = "Ternary plot showing the composition of the variance for each issue. Medicare is the most ideological issue, immigration is the most idiosyncratic, unions is the most unstable.", dpi = 300, fig.width = 6.5, fig.height = 5.5, warning = FALSE}

plot.df <- as.data.frame(overall_decomposition_table)
plot.df$Issue <- rownames(plot.df)
plot.df$Issue <- gsub("_", " ", plot.df$Issue)

### Manual positioning of labels
### Horizontal position
### Default is left-justified
plot.df$hjust <- 1
### Right-justified
plot.df$hjust[plot.df$Issue %in% c("Average",
                                   "Education", "Guns")] <- 0
### Centred
plot.df$hjust[plot.df$Issue %in% c("Contraception", "Unions")] <- 0.5

### Vertical
### Centred
plot.df$vjust <- 0.5

### Down
plot.df$vjust[plot.df$Issue %in% c("Environment", "Unions", "Contraception")] <- 1
### Up
### plot.df$vjust[plot.df$Issue %in% c("Social Security")] <- -0.5

plot.df$fontface <- ifelse(plot.df$Issue == "Average", "italic", "plain")

plot.df$fontcol <- "black"
# plot.df$fontcol[plot.df$Issue %in% c("Health", "Social Security", "Contraception", "Education")] <- "white"
plot.df$Issue <- ifelse(plot.df$hjust > 0,
                        paste0(plot.df$Issue, "  "),
                        paste0("  ", plot.df$Issue))

suppressWarnings(p1 <- ggplot(data = plot.df, aes(Ideology, Idiosyncrasy, Instability), fill = "green") +
    geom_Tisoprop(value=0.5,colour="grey", fill = "green") + 
    geom_Lisoprop(value=0.5,colour="grey", fill = "green") + 
    geom_Risoprop(value=0.5,colour="grey", fill = "green") +
    geom_point(fill = "green") +
    geom_text(aes(label = Issue, hjust = hjust, vjust = vjust, fontface = fontface), size = 3, colour = "black", fill = "green") +
    coord_tern() +
        theme_bw() +
        theme(axis.ticks = element_blank(),
              legend.background = element_blank(),
              legend.key        = element_blank(),
              # panel.background  = element_blank(),
              panel.border      = element_blank(),
              strip.background  = element_blank(),
              plot.background   = element_blank()))
print(p1)

```

```{r,include=FALSE}
# Tabulate low and high knowledge

lowknow_decomposition_table <- decomposition_est_2[,1,]*100
lowknow_decomposition_table <- rbind(lowknow_decomposition_table[sortbyideo,],colMeans(lowknow_decomposition_table)) 
rownames(lowknow_decomposition_table)[nrow(lowknow_decomposition_table)] <- "Average"
colnames(lowknow_decomposition_table) <- c("L:Ideo","L:Idio","L:Inst")
print(round(lowknow_decomposition_table[sortbyideo,]))

highknow_decomposition_table <- decomposition_est_2[,2,]*100
highknow_decomposition_table <- rbind(highknow_decomposition_table[sortbyideo,],colMeans(highknow_decomposition_table)) 
rownames(highknow_decomposition_table)[nrow(highknow_decomposition_table)] <- "Average"
colnames(highknow_decomposition_table) <- c("H:Ideo","H:Idio","H:Inst")
print(round(highknow_decomposition_table[sortbyideo,]))

byknow_decomposition_table <- cbind(lowknow_decomposition_table,highknow_decomposition_table)
byknow_decomposition_table_alt <- byknow_decomposition_table[,c(1,4,2,5,3,6)]

# Test differences between low and high knowledge

decomposition_average_chain_2 <- apply(decomposition_chain_2,c(1,3,4),mean)
decomposition_average_test <- apply(decomposition_average_chain_2[,1,] < decomposition_average_chain_2[,2,],2,mean)
decomposition_diff_int <- apply(100*(decomposition_average_chain_2[,2,] - decomposition_average_chain_2[,1,]),2,quantile,c(0.025,0.975))

```



```{r,results='asis',echo=FALSE, eval = FALSE}
print(xtable(byknow_decomposition_table_alt, digits=0, align=c("l","c","c","c","c","c","c"), caption="Percent of response variation attributable to ideology, idiosyncrasy and instability for each of 13 issues, among low versus high political knowledge respondents."), comment=FALSE, floating=TRUE,table.placement="t",hline.after = 13)
```

```{r,echo = FALSE, fig = TRUE, fig.cap = "Difference in the percentage of opinion variation explained by ideology and idiosyncracy for high-knowledge respondents versus low-knowledge respondents.", fig.width = 7, fig.height = 3, dpi = 300, warning = FALSE, message = FALSE}


plot.df <- as.data.frame(byknow_decomposition_table)
plot.df$Issue <- rownames(plot.df)
plot.df <- melt(plot.df)
plot.df$Respondent <- ifelse(grepl("L:",plot.df$variable),
                             "Low knowledge",
                             "High knowledge")
plot.df$variable <- gsub(".*:","",plot.df$variable)

plot.df$Issue <- gsub("_", " ", plot.df$Issue)
plot.df <- plot.df %>%
    group_by(Issue, variable) %>%
    arrange(desc(Respondent)) %>%
    summarize(Difference = diff(value))

### Arrange the issues so the increase in ideology greatest
issue.levels <- sort(with(subset(plot.df, variable == "Ideo"),
                          by(Difference, list(Issue = Issue), unique)))
issue.levels <- names(issue.levels)
plot.df$Issue <- factor(plot.df$Issue,
                        levels = issue.levels,
                        ordered = TRUE)
plot.df$variable <- recode(plot.df$variable,
                           "Ideo"= "Difference in ideology\n(high knowledge are more ideological)",
                           "Idio"= "Difference in idiosyncrasy",
                           "Inst" = "Difference in instability\n(high knowledge are more stable)")
ggplot(data = subset(plot.df, variable != "Difference in idiosyncrasy"),
              aes(x = Issue, y = Difference)) +
    geom_bar(stat = "identity") +
    coord_flip() + 
    facet_wrap(~variable) +
    theme_minimal()

```

Previous research has suggested that individuals fall into different groups which exhibit different levels of opinion stability [@hill2001extension; @hill2001classification] and opinion structure. Figure 2 shows that, when we split the sample into low and high knowledge individuals, and allow for different ideological dimensions in each half of the sample, we find that high knowledge individuals have more ideologically structured preferences on all but two issues and have less response instability on all issues.  Averaging across issues, ideological structure explains `r round(byknow_decomposition_table[nrow(overall_decomposition_table),4])`% of variation among high knowledge respondents versus `r round(byknow_decomposition_table[nrow(overall_decomposition_table),1])`% among low knowledge respondents, with a 95% posterior interval for the difference running from `r floor(decomposition_diff_int[1,1])`% to `r ceiling(decomposition_diff_int[2,1])`%.  Response instability accounts for `r round(byknow_decomposition_table[nrow(overall_decomposition_table),6])`% of response variation among high knowledge respondents versus `r round(byknow_decomposition_table[nrow(overall_decomposition_table),3])`% among low knowledge respondents, with a 95% posterior interval for the difference running from `r floor(decomposition_diff_int[1,3])`% to `r ceiling(decomposition_diff_int[2,3])`%.  The single ideological dimensions that best predict responses among low and high knowledge respondents are different: among low knowledge respondents, attitudes on guns and unions are the attitudes that most strongly predict other issue preferences, whereas they are among the least predictive among high knowledge respondents. 

# Discussion and Conclusion

These estimates will not come as a great shock to public opinion researchers. 
As we noted at the outset, the ideas that citizens do not organise their views as strongly as
political parties and elites and are sometimes unstable in the preferences 
they express, are hardly novel.  The value of this paper is in providing a 
quantitative decomposition: *how much* of the preference variation we see is due to 
each of these factors?  On average, in a particular survey instrument querying a 
particular set of issues at a particular moment in US history, about $1/7$ of 
cross-sectional variation reflects a common ideological dimension of variation, $3/7$ is real 
but idiosyncratic variation in individuals' views, and $3/7$ is response instability. 

The substantial variation across political issues in the degree to
which citizens have "real preferences" (either ideological or
idiosyncratic), plus the further variation in the extent to which
those preferences are ideological rather than idiosyncratic, indicates
that our results might change substantially with a different sample of
issues.  Our results might also change with different response categories or
different wordings of those response categories -- we note
Broockman's own cautionary remarks regarding the choice of policy
options [@broockman2016approaches, 187] -- but we have no reason to
believe that differently worded response options would
affect our results in any particular direction. The issues chosen include a variety of types of issues (some
primarily economic, some primarily sociocultural), and all of the
issues are sufficiently prominent for researchers to be able to
identify the preferred positions of elite actors
[@broockman2016approaches].  If lower profile issues were chosen
instead, we would expect the estimated ideological and idiosyncratic
shares of the decomposition to decline and instability to increase.

The variation we observe by political knowledge levels suggests that
the results might also vary in a more representative sample. Those who
responded to the second survey wave had somewhat higher levels of
knowledge than those who only responded to the first wave.  This suggests that our findings, if
anything, over-estimate the portion of variation accounted for by
ideology, and under-estimate the portion of variation accounted for by
instability.

In addition to the ways that our estimates are context dependent
because of the structure of the survey data we work with, they are
also model dependent because of the way we structure the analysis.
While the model we specify enables an elegant variance decomposition,
it makes certain functional form assumptions in order to enable that
decomposition.  While we think these are fairly innocuous with respect
to the general conclusions, the exact estimates would surely change
under a different approach to decomposition.  Because the data we work
with has only two waves, we also cannot distinguish instability from
real opinion change, which is a fourth potential source of variation
in response, albeit one that is likely to be more important over
longer time scales than six weeks.  The framework described above can
be easily extended to perform a four-way decomposition given
additional survey waves.  The key to doing so is the logic used by
@converse1964nature: if opinion change is real, opinions in
adjacent waves will be more similar than those in non-adjacent
waves. In the context of our model, this would most naturally take the
form of a normally distributed random-walk process, for each respondent, for each item. We have also set aside
possible opinion processes that would fundamentally undermine the
logic of this sort of decomposition, such as the possibility that
opinion on one issue might predict changes in opinion on other issues
across waves (dynamic constraint).

In sum, our findings put quantitative figures behind the argument made
by @broockman2016approaches that real idiosyncratic preference
variation is much greater than real ideological preference variation
among the US public.  This is important not only for adjudicating
between theories of public opinion, but also because of the widespread
use of unidimensional models, both formal and informal, to understand
interactions betweens elites and the public.  The fact that roughly
3/4 of real public opinion variation does not reflect a unidimensional
spatial dimension, and the fact that a second spatial dimension is
nearly as powerful as the first, does not necessarily undermine the
use of such models.  There is, by construction, always a single
dimension that captures more opinion variation than any other, and
positions on that dimension will therefore always be important to
electoral interactions.  However, conclusions drawn from models
focusing exclusively on a single dimension may or may not be robust to
the fact that this dimension captures only a small proportion of the
real opinion variation on which basis citizens may make choices and to
which elites may have opportunities to respond.


# References

